function phase = phase_canonical_correlation_analysis(data,stimulusFreq,Fs)
%计算SSVEP实验在不同刺激频率下的phase delay
%---input---
%data:size1频率，size2trial，size3数据点，size4通道
%stimulusFreq：刺激频率,本程序仅适配四类问题
%Fs:采样频率
%---output---
%phase:不同频率下的相位差，单位为degree

channel = 15; %Oz

fmin = min(stimulusFreq);
fmax = max(stimulusFreq);
nmin = round(0.08 * fmin - 1);
nmin = (nmin >= 0) * nmin;
nmax = round(0.22 * fmax);
n = nmin:1:nmax;

%类数改变此处需要更改
[n1,n2,n3,n4] = ndgrid(n);
N = [n1(:),n2(:),n3(:),n4(:)];

for i = 1:size(N,1)
    n_list = N(i,:);
    for stimulus = 1:size(data,1)
        for trial = 1:size(data,2)
            x = data(stimulus,trial,:,channel);
            x = reshape(x,size(data,3),1);
            y = fft(x);
            p = unwrap(angle(y));
            L = length(x);
            P2 = abs(Y / L);
            P1 = P2(1:L / 2 + 1);
            SNR((stimulus - 1) * size(data,1) + trial) = narrow_signal_noise_ratio(P1,stimulusFreq(stimulus),Fs);
            deltaf = Fs/size(y,1);
            f_index = round(stimulusFreq(stimulus)/deltaf)+1;
            pe((stimulus - 1) * size(data,1) + trial) = p(f_index) - n_list(stimulus) * 360;
        end
    end
end
